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Abstract 

The perovskitc crystal BaTi0 3 is modeled using a first-principles based effective Hamiltonian and molecular 
dynamics simulations are performed to estimate the pyroclcctric response. The electrocaloric temperature 
change, AT, is calculated for different temperatures and externally applied electric fields. It is found that 
it is possible to achieve a large AT, around 5-6 K, for a relatively small electric field gradient, less than 
100 kV/cm, if the applied fields have a small absolute magnitude. 
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1. Introduction 

A pyroelectric crystal develops a sponta- 
neous electrical polarization when its temperature 
changes. [l[ Its enthalpy can be expressed as H = 
E — £P, where E is the internal energy, £ is 
an applied electric field, and P is the polariza- 
tion. It is possible to cycle the temperature and 
electric field, to drive the crystal through an or- 
der/disorder phase transition, such as the ferro- 
electric/paraelectric transition in perovskite crys- 
tals, to convert between heat and electric energies. 
The phenomena known as the electrocaloric effect 
(ECE), in which electrical energy is used to induce a 
temperature change, is directly related to this ther- 
modynamic cycle. The ECE is not only interesting 
from a fundamental science perspective, but also 
holds great potential for future technologies such 
as solid-state refrigeration. 

The ECE was first discovered in 1930 0, but 
there were few studies pursued until ferroelectric 
crystals, with permanent dipoles, were discovered. 
The initial studies focused on bulk ceramics includ- 
ing Rochelle salt (3, SrTiOa H BaTi0 3 (BTO) [|, 
and doped alkali halides |6j, [D, |8| ; however, the en- 
gineering potential of these crystals was never real- 
ized. 



In 2006 Mischcnko et al. discovered the "giant 
electrocaloric effect" in Zr-rich PbZro.95Tio.05O3 
(PZT) thin films. Q By using a thin film it was pos- 
sible to apply ultra-high electric fields, to achieve an 
ECE AT of 12 K. The discovery of the giant ECE 
has revitalized the field and has inspired numerous 
recent studies investigating the impact of ultra-high 
electric fields on the pyroelectric response of other 
ferroelectric materials. 

The initial work focused on Zr-rich PZT, [<| 
and although this alloy is still being investi- 
gated, [HI [Hi there has also been a great deal 
of interest in other perovskite crystals includ- 
ing (1 - liPb^/aNb^aOa-iPbTiOa 0, E3, 
BTO M,M, M, M 0, BaTi0 3 ~SrTi0 3 
superlattices [20j, and SrBi 2 Ta 2 09 layered ox- 
ides [2l| . In addition there has been much 
interest in the development of polymer based 
ferroelectrics, such as poly(vinylidene fluoridc- 
trifluorocthylene) and poly(vinylidene fluoride- 
trifluorocthylene-chlorofluorocthylene) 



22, 2 



'Corresponding author 
Email address: sbeckmanOiastate.edu (S. P. Beckman) 



In this paper we will focus on the BTO com- 
pound for a variety of reasons: it is a well stud- 
ied archetypical perovskite crystal that is relatively 
easy to produce, it exhibits an excellent pyroelec- 
tric response, it does not contain the toxin lead, 
its ferroelectric properties have the potential to be 
tailored by alloying or the creation of superlattice 
structures. The ECE in BTO has been previously 
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studied as a thin film 
in a multi-layer stack. 



and as a thick film 
Surprisingly, off-the- 



shelf BTO multilayer capacitors have been shown 
to exhibit a reasonable electrocaloric response. [HI 
There also has been a study of the properties of 
BTO nanoparticles with core-shell geometries. [l9| 
Previous theoretical studies of BTO have largely 
relied on thermodynamic models, for example the 
Ginzburg-Landau-Devonshire model. [TBI, [13, [13, 
IS S3, El While these have proven successful, they 
require a large amount of data to accurately param- 
eterize the model. One alternative approach was 
demonstrated by Cao and Li in which they created 
a transverse Ising model (TIM) to study the ECE 
in BTO thin-films. They predict that compressive 
epitaxial strains result in increasing the tempera- 
ture where the ECE is maximum and tensile strains 



reduce the temperature. [14( The TIM model is par- 
ticularly interesting because in addition to the spin 
exchange interactions, taken up to four spins, the 
effect of quantum fluctuations were also included. 
It is reported that increasing the strength of the 
quantum fluctuations results in shifting the ECE 
to higher temperatures. Bellaiche et al. have used 
a molecular dynamics approach to investigate nan- 
odots of Pb (Zro.4Ti .6) O3. They report the time 
dependency of the temperature when an alternat- 
ing electric field is applied. [13] It is discovered that 
it may be possible to tune the ECE in nanodots by 
adjusting the depolarization field, which in practi- 
cal terms would be implemented by modifying the 
surface of the nanodot. 

Following this introduction the theoretical meth- 
ods are presented. The results are presented and 
discussed in section [31 The manuscript concludes 
with a summary in section [4] 



2. Methods 

Here we predict the ECE for BTO using the 
model Hamiltonian, 
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which is derived in Refs. 2J, |25[. The true atomic 
structure has properties that are determined by the 
complex chemical bonding between the atoms, but 
in the model system the complexity is reduced; the 
collective atomic motion is coarse-grained by lo- 
cal softmode vectors u(R) and local acoustic dis- 
placement vectors w(R) of each unit cell at R in 
a simulation supercell as well as the homogeneous 
strain components, (771 , f?2 , ??3 , f]5, %)• The terms 
in the Hamiltonian bear a physical significance: 



Eaa^W and 
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the kinetic energy possessed by the local soft 
modes and the local acoustic displacement vec- 
tors, V scll ({u}) is the local mode self energy, 
V dpl ({u}) is the long-ranged dipole-dipole inter- 
action, V shmt ({u}) is the short-ranged interac- 
tion between local soft modes, ^ elas ' homo (?7 1 , • • •, rj 6 ) 
is the elastic energy from homogeneous strains, 
yeias, mho^^j^ j s ^ e e i as tj c energy from inhomoge- 

nous strains, V coup ' homo ({it}, rji, ■ ■ -,r?e) is the cou- 
pling between the local soft modes and the ho- 
mogeneous strain, V coup,lnh °({u}, {w}) is the cou- 
pling between the soft modes and the inhomogc- 
neous strains, and — Z* J2r £-u(R) is the electric 
enthalpy. Each of the terms in Eq|T] is expanded 
to allow the Hamiltonian to be expressed u sing the 
parameterization given in Refs. [24 . |25[ [26. 27. 28|. 
The same set of parameters in Ref. [28| for BaTi03 
is used in this letter. Values of the 25 parameter 
are also given as Supplementary Data to this let- 
ter. The details of the molecular dynamics method 
are explained in Refs. | 



27J. The method de- 



scribed here is encoded in the software package 
feram, which is distributed freely to the scientific 
community^ 

The simulation uses a supercell of size N = 
L x x L y x L z — 16 x 16 x 16 and the temper- 
ature is scanned from 250 K to 900 K at an in- 
crement of +1 K/step. For each temperature in- 
crement, the system is thermalized for 20,000 time 
steps, after which the properties are averaged for 
80,000 time steps, using a time step of At = 2 fs. 
The simulations begin at low temperature with a 
z-polarizcd initial configuration that is random. As 
described in Ref. (28}, thermal expansion is simu- 
lated by including a temperature-dependent effec- 
tive negative pressure p = — 0.005T GPa. It has 
been shown in Ref. [28j that the parameterization 



lr The MD code "feram" can be downloaded from 
http://loto.sourceforge.net/feram/ as free software. 
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used here offers an improved estimation of ferroelec- 
tric (tetragonal) to paraelectric (cubic) transition 
temperature and c/a ratio as compared to earlier 
models 0, H^j . Therefore, we believe that this set 
of parameters also gives a good estimation of the 
pyroelectric properties. 

The complex interactions, such as the interac- 
tion of domain walls with crystal defects and inter- 
faces that are present in experimental systems are 
not present here. This makes direct comparison be- 
tween computational and experimental results chal- 
lenging; however, these calculations are intended to 
isolate the essential physics of the crystal and the 
pyroelectric response due to the bonds between the 
metal cations and the oxygen. It is anticipated that 
the fundamental behavior observed in these calcu- 
lations will also be present in the experimental sys- 
tems. In principle the predicted ECE response will 
be observed in high purity single crystals of BTO, 
although such experimental systems are challenging 
to prepare. 

3. Results and Discussion 

Simulations are performed for constant electric 
fields ranging from 25 to 310 kV/cm at an incre- 
ment of 5 kV/cm. This temperature range focuses 
on the ferroelectric (tetragonal) to paraelectric (cu- 
bic) transitions, where the largest pyroelectric re- 
sponse should be observed. The z-component of 
polarization, P z , is calculated from the MD simula- 
tion results and is plotted versus T at constant £ z in 
Fig- da). The data are fitted with weighted splines 
to allow for the determination of dP z /dT, which is 
shown in Fig. [ljb) and Fig. [2] Cooling-down sim- 
ulations are also performed and hysteresis in tem- 
perature is found when £ z < 25 kV/cm, as seen 
in the £ z = kV/cm simulations in Fig. 1(a). In 
the simulations, the phase-transition order changes 
from first to second around £ z = 25 kV/cm. To 
focus our theoretical effort on the determination of 
the pyroelectric response without including the ad- 
ditional complexity of the latent heat of transfor- 
mation, only fields £ z > 60 kV/cm are used in the 
present calculations. 

Following Ref. 0] , the electrocaloric refrigeration 
AT is estimated using 

-1 r £ ' 2 BP 

AT =aL T w d£ > < 2 » 

as shown in Fig. [3] For heat capacity, experi- 
mental, temperature and electric-field independent 



value of C v = 5.84 g • cm" 3 x 0.434 J • g" 1 • K" 1 or 
2.53 J • enr 3 • K _1 is used in this estimation (29j ; 
the coarse-graining used in the effective Hamilto- 
nian method makes it difficult to correctly estimate 
the heat capacity. The integral in Eqn. [2] ranges 
over A£ = £2 — £\ . The fields used for this integral 
are specified in the legends of Fig.[3]as, "£1 . . . £2-" 

Because the phase transition becomes increas- 
ingly diffuse as the electric field is increased the 
greatest pyroelectric response, 3P z /dT shown in 
Fig. [1] (b), occurs at lower electric fields. In Fig. [3] 
frame (a) the field gradient A£ is varied and £2 is 
held constant at 310 kV/cm. The magnitude of the 
ECE AT is proportional to the magnitude of the 
field gradient, A£. This is entirely consistent with 
the experimental observations. The temperature 
where the maximum AT occurs is determined by 
the isothermal curve in Fig. [1] that makes the great- 
est contribution within the field range bracketed by 
£\ and £2- For example, over the range £\ = 160 
and £ 2 = 310 kV/cm the T = 550.5 and 600.5 K 
curves have the greatest magnitude in Fig. [2] and 
the AT versus T curve, in Fig. [3] (a), peaks just 
above T = 550 K. 

It is observed by Akcay et al. that not only is 
the field gradient, A£, important but also the ab- 
solute magnitude of the fields £\ and £2- In Fig. [3] 
(b) the field gradient is held constant at A£ = 100 
and £\ is varied from 60 to 210 kV/cm. At larger 
absolute fields the diffuse nature of the phase tran- 
sition results in a decrease in the magnitude of the 
ECE AT as well as an increase in the temperature 
where AT is maximum. As a matter of practice 
this means that the ECE will always occur above 
the Curie temperature. Comparing the results in 
Fig. [3] (a) and (b) it is observed that by keeping £\ 
small it is possible to achieve large AT for small 
A£ , and when £\ becomes large it is necessary to 
have a large A£ to have comparable AT. The im- 
plication of this observation is that compounds with 
antiferroelectric ordering, which require large 
may not be the optimal choice for producing a large 
ECE. 

4. Summary 

In summary, we have demonstrated the use of our 
MD method feram to estimate the ECE as func- 
tions of temperature and external fields. The freely 
available method used here can be parameterized 
using first-principles results alone, and therefore 
can be extended to investigate systems that lack 
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Figure 1: (color online) The calculated pyroelectric response, 
(a) The simulated temperature dependence of the polariza- 
tion along the z-direction at constant external electric field 
is determined. The simulated data is plotted as dots and 
are fitted using weighted splines, shown as lines. £ z = 
kV/cm heating-up and cooling-down results are also plotted 
for reference, (b) The pyroelectric coefficient, dP z /8T, is 
determined from the fitted data. 
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Figure 3: (color online) The predicted clcctrocaloric effect 
estimated from Eqn. [2] The integral range &\ £2 is spec- 
ified in the plot legend, (a) The value of £2 is held fixed 
at 310 kV/cm and £\ is varied to give the A£ the values 
250, 200, 150, 100, and 50 kV/cm. (b) The values of £ x 
and £2 are changed together to maintain a fixed value of 
A£ = 100 kV/cm. 
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Figure 2: (color online) The pyroelectric coefficient, 
dP z /8T, is determined as a function of electric field at con- 
stant temperature from the fitted data in Fig. ^ 
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sufficient experimental data to produce an accurate 
thcrmodynamical model. We have demonstrated 
the effect of varying electric fields on the archetypi- 
cal pyroelectric, BTO, and have found that to pro- 
duce a large ECE for a relatively small field gradi- 
ent, it is preferred to not have an antiferroelectric 
ground state, which will allow the use of modest 
fields. While admittedly the external electric fields 
consider here may be slightly high compared to the 
breakdown voltage of polycrystalline BaTiOa; ad- 
vancements in synthesis techniques, such as vapor 
deposition processes, hold the promise of creating 
affordable high quality ferroelectric thin films that 
can withstand these electric fields. 
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